Mapping pure plastic strains against locally applied stress: Revealing toughening plasticity

The deformation of all materials can be separated into elastic and plastic parts. Measuring the purely plastic component is complex but crucial to fully characterize, understand, and engineer structural materials to “bend, not break.” Our approach has mapped this to answer the long-standing riddle in materials mechanics: The low toughness of body-centered cubic metals, where we advance an experimentally led mitigative theory. At a micromechanically loaded crack, we measured in situ the stress state applied locally on slip systems, and the dislocation content, and then correlatively compared with the occurrence—or not—of toughness-inducing local plasticity. We highlight limitations and potential misinterpretations of commonly used postmortem transmission electron imaging. This should enable better-informed design for beneficial plasticity and strength in crystalline and amorphous solids alike.


Measuring strains: finite strain theory
Deformation gradients F, Fe and Fp may further be separated into strain (stretch & shear) and rotation components, as: Where R is the rotation tensor, U is the right stretch tensor and V is the left stretch tensor. It can be shown (73) that T = 2 and T = 2 , where T denotes the transpose. In the case of infinitesimal deformation, rotation and strain are approximated as summative, as is assumed for the HR-EBSD procedure here: e = + e + , where I is the identity matrix, εe is the elastic strain tensor and ω the lattice rotation. However in the case of plastic deformation, strains may be large, and as such finite strain theory must be applied to deformation measurements by DIC. The Seth-Hill family of deformation tensors describe a variety of methods to evaluate strain (73); the simplest, the Biot strain tensor approximates strain where strains and rotations are small: More useful here, and used for plotting all DIC-based strain data in the present study, is the Green-Lagrangian strain tensor which is similar to simple engineering strain, but further accounts for the second order contributions to strain: Finally, for completeness, true strain may be robustly expressed for large strains using the Hencky strain tensor: Hencky = ln ( ) (5)

Further information on strain mapping methods
In its simplest form, the measurement of elastic strain provides a single, volume averaged value for applied or residual elastic strain, e.g. hole drilling (74), X-ray (75) or neutron (76) diffraction measurement of lattice plane spacing. However, in Materials Science one is often interested in determining local spatial variations in elastic strains and small lattice rotations with high resolution, such as across semiconductor devices (40), and for other crystalline (77) and non-crystalline(11) material classes. Over the past few decades, mostly diffraction-based methods have answered this call: nano-beam synchrotron X-ray diffraction (9), various transmission electron microscopy (TEM) diffraction methods (11,40) and scanning electron microscopy-based high resolution electron backscatter diffraction (SEM HR-EBSD) (10), as well as atomic position imaging in the TEM (78). Accessible regions of interest range from hundreds of micrometres with resolution on the order of 100 nm (79), to mapping elastic strain and lattice rotation about individual dislocations (80,81).
For bulk values, total strain measurement may be achieved by a surface strain gauge or an extensometer, however mapping the local distribution of total strain may currently only be achieved by digital image correlation (DIC) methods -at a surface, by tracking relative movements in a reference pattern (82), and, rarely, by a similar approach in a volume (83). The recent emergence of SEM DIC(13) exploits the increased imaging magnifications of electrons over standard optics, but few studies achieve truly nanoscale strain mapping, nDIC, down to the tens of nanometre resolution (14,15) required to characterise slip lines (37) as this additionally depends on producing a suitably refined nanoscale surface pattern to follow. A known limitation of DIC techniques is that in principle, when applied with a single camera, a 2D projection of surface deformation is calculated; however, recent developments in DIC algorithms (84) demonstrate how this can be overcome in the case of discrete slip line-dominated plasticity, where deformation is constrained to a limited number of known crystal planes.

Loading until the elastic limit
Elastic shear in the plane of the cantilever side surface, εe,12, additionally develops as two extended lobes either side of the notch tip with approximately equal and opposite values. Outof-plane elastic strains (εe,13, εe,23 and εe,33) remain close to the noise threshold at the elastic limit.

Loading beyond the elastic limit
Normal elastic strain along the cantilever axis continues to increase at the notch tip, although εe,22 is higher (more than double, step 6) in the lobe to the left side than that on the right throughout the elastoplastic deformation regime. Contrastingly, in-plane elastic shear strain εe,12 increases more substantially along a ~45° path ahead of the notch tip on the right side, reaching ~1% at step 6, whilst the symmetric feature on the left is limited to ~0.5%. Although the other normal elastic strain components continue to display simple Poisson contraction, and changes to out-of-plane shear εe,23 are negligible, a substantial plume of heightened εe,13 elastic shear develops to the left of the notch; the reason for this is currently unclear.
The slip line operates at an in-plane elastic shear strain εe,12 upwards of 0.5% and the axial total strain εtot,22 generated therein continuously increases as the cantilever is further bent. In addition, a zone of more diffusely distributed total strain, where the strain levels far exceed those achievable by elasticity alone, additionally develops throughout a broad region stretching several microns wide and tall ahead of the notch.

Unloaded state after deformation
After unloading (step 7), all elastic strain components reduce significantly, by approximately 50%, and residual stress is symmetrically distributed about the notch tip; only the out-of-plane εe,13 shear foregoes this trend, retaining the aforementioned heightened positive shear on the left side of the notch only.

Dislocation structures observed by transmission electron imaging
As might be expected, a concentration of linear crystal defects, dislocations, is observed below the notch tip, arranged in a complex network of interlocked tangles. Further, ladder-appearance features of dislocation pile-up extend outwards from the notch tip directed ~45° and ~60° to the horizontal, see skeleton-free TEM image in Fig. S11. The liftout geometry excludes the possibility that these ladders might be any mix of screw to edge type for the four {110} < 111 > slip systems in Fig. 3, whose planes and parallel to the viewing direction and intersect the side surface of the microcantilevers at 45° to the horizontal.

Discrepancy in dislocation density values measured by GND and TEM
The measured ρGND is found to have values approximately double those measured for ρtot; for a comparison on the same colour scale, see Fig. S7. One must note that the HR-EBSD calculation of GND content is sensitive to all dislocation types, including those mostly invisible to the TEM due to edge-on imaging here. However, on the other hand, the GND maps only measure a dislocation content required to fulfil the stepwise lattice rotations, with a 75 nm step here; the undetected, statistically stored dislocations, SSDs, are additionally imaged in transmission. BF-STEM imaging here is useful as contrast is generated from electron loss through (low-angle) coherent scattering by dislocations, whilst (low-to-high angle) mass contrast incoherent scatter and crystal orientation contrast are both minimised as the cantilever is a pure single crystal; hence the most complete dislocation population possible is captured, without significant contribution from other contrast sources. This differs with commonly used weak beam methods for dislocation density analysis, which only capture a specific subset of the total dislocation density; an extinction correction factor is then used (85), which spuriously places assumptions as to the uniformity of the distribution of dislocation types, and hence multiple imaging conditions of a same region may be required.
It is known that there is a minimum ρGND that can be measured, i.e. a noise threshold; this is higher (5.9 × 10 13 m -2 ) with the applied surface Pt speckle pattern here than without for tungsten(25) (5.1 × 10 13 m -2 ). This may be attributed to the base-level noise in the measurement of ω rotation components, leading to a non-zero Nye tensor even before loading, which worsens with further mapping due to accumulating surface contamination. Similarly, the TEM intercept method is limited at the upper end of the range, when the dislocation array is highly dense; it is difficult to capture closely spaced dislocations in the skeletonisation. This is seen in the vicinity of the notch tip, Fig. 4C, where the measured ρtot is maximised at around 2 × 10 15 m -2 .
Despite these differential lower and upper bounds, ρGND is nevertheless higher than ρtot throughout much of the conventional notch-tip plastic zone. This may be because ρtot reduces upon TEM lamella fabrication, by dislocation migration to the surfaces. However, unlike the low lattice resistance of some metals such as copper which results in room temperature annealing in short timescales (86), the high Peierls stress of tungsten (87) is thought to resist this annealing process. Alternately, one may also consider that the cantilever bending geometry here inherently leads to elastically-mediated lattice rotation preceding any local dislocation activity; hence a phantom elastic component may additionally contribute to the measured Nye tensor, artificially further raising the inferred GND density, as previously seen (25). Nevertheless, upon unloading ( Fig. S7) this background contribution should mostly disappear. Finally, one may note in Fig. S7 that for the low dislocation density, elastically relaxed region of the cantilever root, far to the left of the notch, the logical ranking of ρtot > ρGND is adhered to, despite this being at the maximum accumulation of surface contamination from repeated scanning. This observation suggests that the storage of dislocations itself is a major contributor, causing the Kikuchi pattern quality to diminish, and in turn the noise-induced fictitious GND density to increase. This is a universal concern for the application of HR-EBSD. However, it does not contradict key observations here, such as the symmetry of ρGND about the notch.
6. Note on the apparent asymmetry of the notches One should note that the apparent asymmetry of the notch shape in elastic and total strain maps, Fig

Extraction of the pure plastic component of deformation, Fp
In theory, an obvious numerical approach to extracting the pure plastic deformation gradient tensor from the DIC and HR-EBSD data is to inverse the elastorotational gradient Fe, and use this to left-multiply the total deformation gradient: In fact, it is expected that in most cases for engineering materials, the plastic strains that must be sustained for these to be considered 'damage tolerant' are sufficiently large that the contribution of elastic strains to the total are negligible(88); however, the lattice rotations are not.
In practice, one must consider the configurations in which the data is collected. For example, if one executes the digital image correlation procedure in the forward sense, that is the image of the deformed condition relative to that of the initial undeformed case (i.e. calculation of F), the output may be considered on either the regular grid of the initial image, or on the irregular, deformed grid, as illustrated in Fig. A.1 of (8). In Fig. 2 of the present study, the total strain data from F is plotted on the deformed grid, for example. For the extraction of Fp, it is convenient if both the DIC and HR-EBSD data are plotted on a regular grid, with ideally the same pitch, to facilitate overlay.
From Fig. 1, it is apparent that the configuration common to both the total and elastorotational deformation is B, the final configuration. Therefore, this was chosen as the configuration on which to perform the calculation of Fp, noting that Fp would therefore be plotted on a configuration corresponding to neither of its 'natural' states -i.e. neither B0 nor Bp. In theory, it is possible to return the regular grid of Fp data to the B0 or Bp configurations through a reverse transformation of the gridpoint positions using F -1 or Fe -1 , respectively. This was not considered necessary in the context of the present study.
In order to generate a regular grid of datapoints in configuration B whilst avoiding introducing data processing artefacts and compounding error (i.e. limiting unnecessary data interpolation) the reverse digital image correlation between image pairs (initial and step X) was performed. Hence F -1 was instead calculated, based on a regular grid of correlation subsets in the deformed configuration at loading step X. Similarly HR-EBSD scans produced, by default, a regular grid of datapoints in the deformed configuration. With this in mind, it is more judicious to consider an alternative expression of Eq. 6a, as used in the present study for the extraction of Fp: as this minimises the number of inversion operations required compared to Eq. 6b, which would require the inversion of both measured F -1 and Fe independently. Fe was compiled according to: where I is the identity matrix, εe is the elastic strain tensor and ω the lattice rotation, as in Fig. 2. F -1 was compiled as: where δij is the Kronecker delta and u' is the displacement field that relates the position of elements in the initial state, relative to their position in the deformed state. dx is a line element in the deformed state, and dX (Fig. 5) equivalently in the initial state, i.e. dX = F -1 dx, as in (8).
As the SEM-DIC strain mapping method only allows the in-plane strain components to be calculated, Fp was similarly restricted: therefore only the four 1,2 indices components of Fe were used for the calculation of Fp.
The magnifications for SEM imaging of the Pt surface speckle pattern were chosen such that one of the common binary subset sizes would be 75 nm wide, or simple fractions or multiples thereof, e.g. for the smaller cantilever size, 16 px = 75 nm. Where necessary regridding of the F -1 field could then be performed using bi-cubic interpolation based on standard Matlab functions.
9. An annealed slip band The concept of an annealed slip band, along which mobile dislocations move rapidly before exiting the sample at a surface is supported by the absence of reverse plasticity along this line upon unloading of the cantilever and the low total dislocation density in this location. Unlike the distributed dislocation array captured by GND mapping, most dense at step 6, Fig. 4, but of which some nevertheless remains upon unloading, the very large operative dislocation population responsible for the intense right-hand slip line has been lost -reverse plasticity there would require significant further nucleation of dislocations, for which stress accumulation upon unloading is apparently insufficient to motivate it. Indeed, one should consider that the, at least, 700-strong population of mobile dislocations identified in Fig. S10 to generate the southeastern plasticity could not simply be stored within the cantilever, as this would equate to an average dislocation density of 2 × 10 15 m -2 along the entire slip band, which is clearly not measured by TEM, Fig. S11. 10. Analysis of results based on finite element crystal plasticity theory, and other models Within finite element theory, that is the step-wise partitioning of space, as applied to crystals(73) (e.g. CPFEM), dislocation populations may be subdivided into geometrically necessary dislocations (GNDs) and statistically stored dislocations (SSDs, density ρSSD). Whilst the former fulfil the geometrical incompatibilities in the plastic deformation of adjacent cells, i.e.
Curl F p ≠ 0, measured as rotation of the crystal lattice in the step-wise HR-EBSD mapping, the latter are arranged in such a fashion that no net lattice rotation is generated between neigh-bouring cells by this population. This classification of dislocations is not fundamental; however, it does serve in the understanding of the observations here of GND and total dislocation contents, and how this may relate to plasticity. The fact that after unloading the 45° band of intense plasticity and high ρGND on the right side of the notch tip displays no salient features on the ρtot map, whereas the low-plasticity, high ρGND 45° band on the left presents a raised total dislocation density, indicates a disparate distribution of SSDs. Hence ρSSD in the lefthand band must be higher than on the right side of the notch tip. Accurate quantification of this ρSSD difference is prohibited for reasons in Supplementary Material 5. It is interesting to note that within this finite element formulation, the population of mobile dislocations m is proportional to the density of immobile, statistically stored, dislocations SSD for a given slip system α (via a slip system-dependent interaction parameter) -the contribution of GND content to m remains spatially symmetrical about the notch throughout testing. At the same time, the net shear rate ̇= m and therefore depends on the dislocation velocity v α which for the present discussion varies of significance as ∝ sinh� � − � P + m ��, hence: where for positive applied shear stress τ α on the slip plane, µ is the shear modulus, b the Burgers vector of slip, P is the density of dislocation components parallel to α of all slip systems, which will hence induce a passing resistance to the movement of the mobile α dislocations that constitute m , and A and B are positive constants for fixed temperature and slip system type (89). Hence the effect of m on the strain rate is not simply determined as decreases with increasing m .
The present single crystal example demonstrates a certain non-physicality of this finite formulation: although a low, fast moving mobile dislocation population is experimentally found to be beneficial to plasticity, this low m serves to repress the programmed multiplication processes in the CPFE formulation, e.g. dipole formation (89), whilst no dislocation loss at free surfaces or other sinks is usually envisaged either. The experimental result here would require a highly active population of dislocation sources along the dominant slip line in the CPFEM context -in reality, dislocations are likely mostly generated in the high stressed proximity of the notch tip. Comparison relative to existing CPFEM studies (34,42) of microcantilever bending using non-local formulations that allow GNDs and hence a length-scale to be introduced yields discrepancies concerning the role of GNDs. The storage of GNDs in such simulations (34,42) has been interpreted to impose an additional resistance to the motion of the mobile dislocations, resulting in small-scale strengthening effects. However, the present GND density symmetry and plastic asymmetry would suggest this description is incomplete. In such CPFEM cases, the nucleation rate and spatial distribution of sources, as well as the possible starvation of mobile dislocations, are not considered.
Discrete dislocation models, which do account for dislocation nucleation and loss of the individual dislocations, have similarly been applied to microcantilever bending (28); equivalently to GNDs, nucleation source density was found to strongly control small scale strengthening effects. However, the simulated Frank-Read sources operated when a threshold stress was locally exceeded for a prescribed time; hence a variable rate of activation of the source has not been robustly considered.
11. Discussion points concerning loading-or strain-rate It was necessary, here, to pause in situ loading (otherwise at 0.05 MPa m 1/2 s -1 for the 5 µm height cantilever, 0.06 MPa m 1/2 s -1 for the 3 µm one; these are relatively slow rates (22)) for 20 min at each strain mapping step for diffraction map and speckle image acquisition, under load: effectively a variable loading rate fracture test. In theory, for BCC tungsten at room temperature, some dislocation motion can occur during these holding periods; nevertheless, the coherency in J-integral against crack length when comparing with continuous loading (Fig. S1) indicates the pauses have little effect on the deformation mechanistics and cracking during the dynamic phases.
12. Discussion points concerning the toughness of tungsten Current understanding of fracture processes across brittle to ductile transitions (BDT) find consensus in crack tip plasticity in crystals being a two event process (90). For effective toughening, it requires both efficient dislocation nucleation and rapid glide i.e. high dislocation mobility within the plastic zone (44). Previous seminal work on tungsten (45) found nucleation to be limiting only well below the BDT, by cryogenic testing after different amounts of cold work-induced dislocation source arrays were introduced. In the transition regime, and at higher temperatures, dislocation mobility controlled crack tip toughening processes.

Size effects study
Smaller notched cantilevers, 3 µm in height, were equivalently deformed, Fig. S9. The conditional fracture toughness was lower for these reduced dimensions, 5.1 MPa m 1/2 , as was the maximal angle of rotation of the beam relative to the root before rapid crack growth. The magnitude of plastic strain achieved below the notch was similarly reduced compared with the larger cantilevers; however, the asymmetry of the distribution of strains and lattice rotations relative to the notch was maintained, albeit more subtly. Plasticity is observed to focus on the right side of the notch, whilst the normal stress along the cantilever axis is a maximum ahead of the notch on the cantilever root side. Finally, the region of heightened GND density forms a vertical band, following the crack propagation direction, and displays little lateral spread, unlike the notch plume in the 5 µm cantilever.
Although elastic stresses and strains are similarly distributed during loading of both notched cantilever sizes, the lateral distribution of geometrically necessary dislocations is considerably lesser in the smaller sample and achieves lower maximum densities, yet remains symmetrical about the notch: dislocation multiplication initiates at the notch tip and GNDs extend only as far as soft pile-up at the axial tension-compression boundary (see σ22 component) permits. Fig. S1: Comparison of strain mapped and continuous deformation cantilevers.

Supplementary Figures
Comparison of current notched cantilever loading here (larger cantilever size, Fig. 2, blue) with previous data (grey) acquired without interruptions, showing good agreement in both (A) loading curves and (B) J-integral crack growth lines. In (A), the mode I stress intensity factor at the crack tip calculated based on linear elastic fracture mechanics (LEFM) (22) is plotted against displacement of the sphero-conical tip contact point. Note that although this LEFM assumption is only truly valid until the elastic limit (step 2), it serves as an effective method to compare crack tip stress states, normalising any differences in exact cantilever geometry between the two testpieces, such as in loading arm length. (B) confirms that no measurable crack growth occurs during the displacement holding periods for acquisition of EBSD maps and speckle pattern images. The conditional fracture toughness measured using the plasticityconsidering J-integral method for the larger cantilever size here, 6.8 MPa m 1/2 at 0.05 MPa m 1/2 s -1 , lies comfortably within the 6.6 -7.0 MPa m 1/2 range previously determined (24) for stress intensity factor rates of 0.03 -0.1 MPa m 1/2 s -1 .  Table S1 for both cantilever sizes. The scale bar in (A) is 50 µm long; the scale bars in (B) are 5 µm long.   (A) Atomic force microscopy scan of the root-notch region of a large size cantilever, equivalent to that in Fig. 2; it was not possible to scan the exact same cantilever due to the eventual damage caused by AFM scanning, which should not compromise the possibility of obtaining the liftout sample for TEM imaging, Fig. 5. The scan step size was 20 nm. The Pt surface speckle pattern is visible across the entire scan area. Profiles (black and orange) averaged across a 20 px width are extracted on the right; the height of the step across the slip line (profile 1) was measured to be 7.5 nm. A second degree polynomial was fit to the profile segments as a guide to the eye; non-planarity, particularly profile 2, was due to the unavoidable FIB milling artefact near the notch. The white patch mid-right is surface contamination. (B) For the larger size cantilever in (A), the difference between the pixel values of the initial (undeformed) and final (unloaded, after deformation) SEM images of the sample surface is plotted. The slip line extending at ~45° south-east from the initial notch-tip position is clear; no equivalent line is observed here, or by AFM (profile 2), extending at ~45° south-west of the notch tip. A 5 × 5 px 2 Gaussian smoothing kernel was applied to the differential image. The dashed white rectangle in (B) delimits the AFM scanned area in (A).
(C) Total normal strain map along the cantilever axis for the same cantilever as in (A,B); the final (eighth), unloaded step is presented here. The slip line heading southeast from the notch tip is again visible. The map is lower resolution than those in Fig. 2 due to increased imaging noise: the DIC subset size was 169 × 169 nm 2 (50% overlap). The scale bars are 1 µm long.  Considerably lower total dislocation densities are measured below the notch than GNDs, indicating either an extensive loss of dislocations upon thin film extraction, or an overestimation of GNDs in the energy-minimisation of the Nye tensor solution method (30,31). The scale bars are 1 µm long. and confirms the direction of the slip step (see Fig. S5), whilst also suggesting opportunities for future improvements of the correlative strain mapping setup, such as stage tilt back-andforth at each loading step to image the speckle pattern normal to the cantilever surface of interest. The pitch of the square grid for Fp maps is 75 nm. All scale bars are 1 µm long.  Fig. 2, the axial normal stress along the cantilever length and the maximum total shear strain is measured at several steps during loading. Between loading steps 3 and 4, rapid crack growth occurs; relatively little plastic deformation is observed to accumulate prior (strain scale retained from Fig. 2 for comparison), although asymmetry of plasticity is again observed at this size scale, being more intense southeast of the notch, arrowed in red. (B) The accumulation of geometrically necessary dislocations ahead of the notch is similarly measured to be more conservative than the larger cantilever, both in lateral spread and intensity. This coheres well with the lower measured toughness for the 3 µm cantilevers, resulting from lesser crack tip plasticity. The spatial resolution for HR-EBSDbased data is 75 nm; total deformation was presently calculated with 24 × 24 px 2 subsets with 25% overlap, yielding an 84.4 nm subset spacing. Errors of measurement for all quantities are again given in full in Methods. All scale bars are 1 µm long.  strain mapped surface was imaged by bright field scanning TEM. This aimed to capture the most dislocations possible (in contrast to alternative weak-beam methods that selectively exclude certain dislocation orientations). A higher resolution image of the notch tip region is given in (C). Diffuse background contrast variations result from lifting the lamella close to the FIBbed surface, rather than unoptimised preparation for TEM. Following manual skeltonisation of the dislocation array (B,D), two masks for intercept measurement were created (E) with horizontal and vertical lines, or 45° diagonal lines. The mask size was 21 × 21 px 2 = 75 × 75 nm 2 , i.e. equal to the HR-EBSD scan step. Intercepts between the masks and the dislocation skeleton were calculated using a Matlab code; from this the dislcoation density maps (F) were produced for both masks. The close similarity of the two ρtot maps demonstrates the robustness of the method relative to the choice of mask. The scale bars are 1 µm long; that for the diffraction pattern is 4 nm -1 .

Fig. S12: Error values for the Fp components.
Based on the standard noises of measurement determined for HR-EBSD and DIC data, the propagation of uncertainty into the calculation of Fp was performed according to (91). The matrix inversion procedure results in higher absolute error along the line of extensive slip than in the remainder of the cantilever; the relative error is lower however along the slip line for  Smaller, Suppl. Fig. 9 3.3 2.9 1.1 9.4